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Abstract 

Genetic evaluation using animal models or pedigree-based models generally 
assume only autosomal inheritance. Bayesian animal models provide a flexible 
framework for genetic evaluation, and we show how the model readily can 
accommodate situations where the trait of interest is influenced by both auto- 
somal and sex-linked inheritance. This allows for simultaneous calculation of 
autosomal and sex-chromosomal additive genetic effects. Inferences were per- 
formed using integrated nested Laplace approximations (INLA), a nonsam- 
pling-based Bayesian inference methodology. We provide a detailed description 
of how to calculate the inverse of the X- or Z-chromosomal additive genetic 
relationship matrix, needed for inference. The case study of eumelanic spot 
diameter in a Swiss barn owl (Tyto alba) population shows that this trait is 
substantially influenced by variation in genes on the Z-chromosome 
{a\ = 0.2719 and a\ = 0.4405). Further, a simulation study for this study sys- 
tem shows that the animal model accounting for both autosomal and sex-chro- 
mosome-linked inheritance is identifiable, that is, the two effects can be 
distinguished, and provides accurate inference on the variance components. 
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Introduction 

In general, quantitative genetic methods implicitly assume 
only autosomal inheritance when estimating variance 



components and heritability for different types of traits 
(Qvarnstrom et al. 2006; Foerster et al. 2007; Forstmeier 
et al. 2011). In this study, we explore the consequences of 
not modeling sex-linked inheritance when estimating 
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additive genetic effects when some of the genes control- 
ling a trait are located on a sex chromosome. 

The heterogametic parent, for example, XY males in 
mammals, and ZW females in birds, only gives its X/Z 
sex chromosome to its homogametic offspring (i.e., XX 
females in mammals and ZZ males in birds). Hence, 
when the selection acts strongest on the heterogametic sex 
in the population, the genes on the X/Z sex chromosome 
will be exposed to selection only half of the time com- 
pared with genes on autosomes (Rice 1984). Thus, selec- 
tion will influence genes located on the autosomes the 
most (Charlesworth et al. 1987), and as a result, we 
would expect to see a much slower change over time in 
genes on the sex chromosome than in genes located on 
the autosomes. How natural selection and sexual selection 
affect the evolution of a trait will depend on whether the 
contributing genes are on autosomes or sex chromo- 
somes. The importance of determining whether a given 
gene, quantitative trait locus (QTL), or part of genetic 
variation that contributes to phenotypic variation is 
located on autosomes or sex chromosomes has been 
emphasized in many studies (Charlesworth et al. 1987; 
Mank and Ellegren 2009; Blackburn et al. 2010). This 
knowledge is especially important in understanding the 
evolution of sexual dimorphism (Rice 1984), but may also 
affect the rate and direction of phenotypic evolution in 
general (Lande 1980; Kirkpatrick and Hall 2004). 

A generalized linear mixed model that offers a powerful 
approach to estimate genetic variance components, such 
as autosomal and sex-linked additive genetic variance, is 
the so-called animal model (Kruuk 2004). In contrast to 
simpler methods such as parent-offspring regression or 
sib designs, animal models utilize information from dif- 
ferent relationships between individuals in large and com- 
plex pedigrees simultaneously. Animal models express the 
phenotypic value of a given trait as a linear sum of fixed 
and random effects, where the different random effects 
have a specified covariance structure. The most important 
structured random effect is the additive genetic effect 
(breeding value), which has a covariance structure given 
by the additive relationship matrix (Lynch and Walsh 
1998). Including the additive genetic effect allows for esti- 
mation of important genetic parameters such as additive 
genetic variance and heritability. 

However, the covariance structure of the breeding val- 
ues reflects a mode where the genetic relationship 
between relatives of the same degree is assumed equal 
irrespective of sex, and as such it corresponds to an auto- 
somal mode of inheritance (in that each individual inher- 
its one half of its autosomal genes from each of its 
parents). This representation of the additive genetic effect 
does not take into account that sex-chromosomal genes 
might contribute substantially to the total additive genetic 
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effect and variance for certain traits. For example, sex- 
linked effects are found in Drosophila (Cowley et al. 
1986), in birds (Sastre et al. 2003), and humans (Pan 
et al. 2007). 

The assumption of only autosomal inheritance may not 
only prevent one from gaining important knowledge 
about where the genes contributing to phenotypic varia- 
tion are located, but may also result in inflated estimates 
of what should be interpreted as autosomal additive 
genetic variance. The latter occurs due to the similarities 
in the inheritance patterns (Grossman and Eisen 1989; 
Lynch and Walsh 1998), and thus the covariance struc- 
ture of autosomal and sex-linked genes. Erroneously 
assuming that all additive genetic variance is due to genes 
on autosomes may result in biased predictions for the 
rate and direction of adaptive evolution (Lande 1980; 
Kirkpatrick and Hall 2004). 

To separate autosomal from sex-linked additive genetic 
variances using the animal model, we need to explicitly 
model sex-linked effects by utilizing the corresponding 
covariance structure of genes on the sex chromosomes. 
The theory on how to construct the necessary covariance 
matrix for inclusion of Z-linked additive effects is pre- 
sented in Fernando and Grossman (1990). 

However, only a few authors within evolutionary quan- 
titative genetics have considered sex-linked additive 
genetic effects within the animal model framework. Fairb- 
airn and Roff (2006) suggested to use the animal model 
for estimating genetic variance due to sex-linked genes in 
the context of evaluating of sexually dimorphic traits, yet 
they did not present any results from the proposed 
model. An extensive version of the animal model was pre- 
sented in Meyer (2008), which, among other genetic and 
environmental effects, also accounted for sex-linked addi- 
tive effects. They used simulated data on an experimental 
design to estimate the variance components using 
restricted maximum likelihood (REML) methods, and 
their results showed that the model was able to disentan- 
gle additive genetic variances caused by sex-linked and 
autosomal effects. To the best of our knowledge, Roulin 
et al. (2010) and Husby et al. (2013) are the only authors 
who have applied an animal model accounting for both 
autosomal and sex-linked additive effects to empirical 
data from natural populations. Roulin et al. (2010) esti- 
mated autosomal and sex-linked heritabilities of a mela- 
nin-based plumage trait (i.e., the size of black spots 
located of the tip of feather of the ventral body side) in a 
wild population of Swiss barn owls (Tyto alba) and found 
that this trait was significantly influenced by sex-linked 
genes. Husby et al. (2013) estimated autosomal and sex- 
linked heritabilities (and additive genetic variances) of 
both morphological and (assumed) sexually selected traits 
for comparison in two long-term (pedigree) studies of a 
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natural population of collared flycatchers (Ficedula albi- 
collis) and a captive population of zebra finches (Taenio- 
pygia guttata). Most traits in both species were not 
significantly influenced by sex-linked genes or showed 
low levels of sex-linked genetic variation. However, wing 
patch size in collared flycather (known to be under sexual 
selection) showed a higher level of sex-linked genetic vari- 
ation. 

The main focus in this article is to show how to 
explicitly model the additive effect of genes residing on 
the larger sex chromosome, that is, the X-chromosome 
which is found, for example, in most mammals and 
some insects (e.g., Drosophila) and the Z-chromosome 
found in birds, butterflies, moths, and some fishes 
(Russel 2006). 

In this study, a simulation study is conducted for the 
barn owl study system to assess the identifiable properties 
of the model assuming both autosomal and sex-linked 
effects and to evaluate the consequences of using a model 
which only assumes autosomal inheritance when the trait 
under study is actually influenced by sex-chromosomal 
genes. We also present a detailed description on how to 
obtain the relevant precision matrices (inverse covariance 
matrices) required to explicitly account for and model 
sex-linked additive effects and set up an extended animal 
model. The objective is to provide a consistent framework 
allowing for estimation of both autosomal and sex-linked 
additive genetic effects using an animal model. 

The methodology presented is also illustrated by 
analyzing the same melanin-based trait as in Roulin 
et al. (2010). Our approach do, however, avoid the 
numerical problems in inverting the precision matrix 
accounting for sex-linked additive genetic effects that 
were reported by Roulin et al. (2010), resulting in more 
precise estimates. 

All inferences in this study are carried out using 
Bayesian methods. One of the main advantages of Bayes- 
ian methods compared with the more traditional REML 
methods is the more accurate representation of uncer- 
tainty in parameter and random variables estimates. 
Bayesian methods allow uncertainty to propagate 
through the model such that all available information is 
contained in the posterior distribution of the parameter 
and random variables in question. Although well estab- 
lished in the field of animal breeding (e.g., Sorensen and 
Gianola 2002), the use of Bayesian methods to tackle 
evolutionary questions has only recently been introduced 
(Kruuk et al. 2008; O'Hara et al. 2008; Ovaskainen et al. 
2008; Hadfield 2010; Steinsland and Jensen 2010; Holand 
et al 2013). We follow Holand et al. (2013), and use the 
Bayesian approximation methodology integrated nested 
Laplace approximations (INLA) introduced by Rue et al. 
(2009). 
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Materials and Methods 
Field data 

We use field data from a wild population of Swiss barn 
owls, a medium-sized nocturnal bird, in western Switzer- 
land. In this study area covering 190 km 2 , 20-80 pairs of 
barn owls breed each year in 110 nest boxes put up in 
barns. We consider the plumage trait diameter of black 
eumelanic spots found on the tip of feathers on the owls' 
ventral side. The data in our study were recorded in the 
period between 1996 and 2007. 

The number and size of eumelanic spots varies both 
within and among populations, and also within families 
(Roulin 2004). Spot size a sexually dimorphic trait, 
where females display on average larger black spots than 
males (females; mean = 13.13 mm [SD = 3.40], males; 
mean = 9.36 mm [SD = 3.96]). In this Swiss population, 
spot diameter has been shown to harbor a high herita- 
bility {h 2 = 0.61) with some variation explained by 
genes on the Z-chromosome (Roulin et al. 2010). Fur- 
thermore, it has been shown that females, but not 
males, are positively selected for large spots (Roulin 
et al. 2010). 

As extra-pair paternity is rare in the barn owls (Roulin 
et al. 2004), a pedigree was constructed by assuming that 
the social parents are the biological parents. Sex of nes- 
tlings was found using sex-specific molecular markers 
typed in blood cell DNA, and from the presence of a 
brood patch in breeding females (Roulin et al. 1999). 

The pedigree consists of AL = 2999 barn owls, with 
1550 females and 1449 males. Plumage spots are 
expressed already at the nestling stage, and spot diameter 
is measured for most individuals in the pedigree 
(N d = 2543, 1333 females and 1210 males). The spot 
diameter data are standardized to have mean 0 and vari- 
ance 1. Further, sex and hatch year is available for all 
individuals in the pedigree and has been found to be 
important for both variation in and selection on plumage 
spot diameter (Roulin et al. 2010). The plumage spot 
diameter is approximately Gaussian distributed (see Fig. 
SI). Mean spot diameter for both sexes and for females 
and males separately for each cohort (i.e., hatch year) is 
given in Fig. S2 and suggests changes in spot diameter 
over the study period. For a more thorough description 
of the fieldwork and methods, study area, and genetic 
analyses, see, for example, Roulin et al (2010) and refer- 
ences therein. 

There are some differences in the dataset used in this 
study compared with Roulin et al. (2010), and the reason 
for using slightly different datasets is farther explained in 
the Discussion. First, Roulin et al. (2010) used a pedigree 
consisting of N p = 3264 individuals captured between 
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1987 and 2007, and this is the same pedigree as in our 
study except that the parents of the individuals with 
hatch year 1996 was included (i.e., they were the "foun- 
ders") in their study. These individuals were also included 
in the dataset used in Roulin et al. (2010), giving 
= 2711 individuals. Second, the dataset in Roulin 
et al. (2010) included all individuals alive in 2007, 
whereas in the dataset used in our study all individuals in 
the pedigree with hatch year before 2006 were taken out. 
Furthermore, in this study, the phenotypic trait was stan- 
dardized to have mean zero and variance 1. In contrast, 
Roulin et al. (2010) standardized the data within each 
sex, which resulted in phenotypic data with mean zero 
and variance 1 within each sex. 

Animal model 

To introduce the animal models, we review models pre- 
sented in, for example, Lynch and Walsh (1998), Sorensen 
and Gianola (2002) and Kruuk (2004). We first introduce 
a Gaussian version of the animal model for only autoso- 
mal loci. The model is further extended to also include 
sex-linked inheritance. 

An animal model for autosomal inheritance (Al) is a 
(generalized) linear mixed model where the observed trait 
values y b i = I,.. .,N d are given by: 

yi = /? 0 + zjp + a ; + e;, (1) 

where /3 0 is an intercept, /? = ft N ) are referred to 

as fixed effects that account for group-specific effects such 
as, for example, sex and hatch year (although in theory 
all Bayesian parameters are random) and zf is a known 
incidence vector. The a/s are individual additive genetic 
effects and are genetically linked random effects also 
known as breeding values. e ; is individual i's residual 
effect, and is an unstructured Gaussian random effect, 
often called the environmental effect in quantitative 
genetics. The parameters /}, e, and a are assigned inde- 
pendent Gaussian priors, /?~ Af(0, t^J), the residual 
effects e~ Af(0, o 2 I), where J is the identity matrix and 
o\ is generally referred to as the environmental variance. 
The additive genetic effects of the autosomal loci are for 
the population, a = (a l ,a 2 ,...., a Np ), assumed to have a 
covariance matrix Aa 2 a , with a dependency structure cor- 
responding to the pedigree u~jV(0, o\A), where A is the 
relationship matrix whose elements are twice the coeffi- 
cient of co-ancestries between relatives for autosomal loci, 
and a 2 a is the additive genetic variance in the base popula- 
tion (see e.g., Lynch and Walsh 1998; Sorensen and Gian- 
ola 2002). According to A, an individual receives half of 
its autosomal genes from each of its parents irrespective 
of sex (Quaas 1976), and a 2 is an estimate for additive 
genetic variance for autosomal loci. Hence, the model in 
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eqn (1) models the additive effects of genes located on 
autosomes. 

To include the additive genetic effects of the sex chro- 
mosomes, we model the additive genetic effect of genes 
residing on the largest of the sex chromosomes, for birds 
the Z-chromosome, and assume the smallest chromosome 
(here W) is inert with respect to additive effects 
(Fernando and Grossman 1990; Ellegren 2007). The total 
additive genetic effect is then partitioned into the sum of 
additive effects due to autosomal genes and additive 
effects due to Z-linked genes. Statistically, it is straightfor- 
ward to include a new random variable in the animal 
model, such as the Z-linked additive genetic effect, given 
its corresponding covariance structure is available. We 
can extend the Al animal model in eqn (1) to an autoso- 
mal and Z-linked inheritance (AZI) animal model 
accounting for both autosomal and sex-linked additive 
genetic effects: 

y { = P 0 + zjp + en + z, + E h (2) 

where z ; is the individual i's additive genetic effects due 
to genes on the sex chromosome. The additive genetic 
effects of the Z-chromosome for the population 
z = {z\,Z2, .... } Ztf f ) are assumed to have a covariance 
matrix Zo 2 z , with a dependency structure corresponding 
to the pedigree and the sex of individuals in the pedigree. 
It is given a Gaussian prior z~ jV(0, g\Z), where Z is a 
matrix whose elements are functions of the coefficient of 
co-ancestries between relatives for the Z-chromosomal 
loci, and a 2 is the variance of additive genetic effects for 
sex-chromosomal genes for the homogametic sex, here 
males, in the base population (Fernando and Grossman 
1990). 

The underlying theory for computation of Z rests on 
some assumptions. The population is assumed to be in 
gametic equilibrium, the additive genetic effect for the 
same allele is assumed to be equal for males and females 
(no dosage compensation Ellegren et al. 2007b; Itoh et al. 
2007), and allelic frequencies are equal in the two sexes. 
Z differs from A because the sex-linked genes are trans- 
mitted in a different pattern than the autosomal genes. A 
female (the heterogametic sex, here ZW) receives all of 
her Z-linked genes z from her paternal parent (the homo- 
gametic parent, here ZZ) and no Z-linked genes from the 
maternal parent (the heterogametic parent, here ZW), as 
mothers pass on their W-chromosome to daughters. On 
the other hand, a male will receive z m from his maternal 
parent and z p from his paternal parent, z = z m + z„. Thus, 
the additive Z-linked genetic variance for noninbred 
males (homogametic sex) is a 2 zm = Var(z m + Zp) = 
Var(z m ) + Var(zp) = a 2 , while for noninbred females 
(heterogametic sex), a \f = Var{z p ) = {l/2)a 2 z . Hence, 
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for a noninbred population with an 1:1 sex ratio, the total 
variance in the population due to Z-linked inheritance is 

< pop = (3/4K 2 . 

Throughout, A will be referred to as the autosomal 
relationship matrix, and Z will be referred to as the 
Z-linked relationship matrix. We assigned inverse Gamma 
priors a 2 * ~ invGamma(a*, b*), where a* = 1 and 
b* = 0.001 to variance parameters <t^, a 2 a , a 2 , and a 2 . 

Modeling Z-linked inheritance - INLA and 
computational issues 

We use integrated nested Laplace approximations (INLA) 
to estimate variances (o^, a^, o^, o^), individual breeding 
values (flj, Z;) and DIC from Al animal models (eqn 1) 
and AZI animal models (eqn 2). INLA is a fast and deter- 
ministic nonsampling-based approach to Bayesian infer- 
ence available for latent Gaussian Markov random field 
(GMRF) models (Rue et al. 2009). It has been shown that 
the Al animal model falls within the class of GMRF mod- 
els, and INLA can be used as inference method (Steins- 
land and Jensen 2010; Holand et al. 2013). 

For INLA methodology to work efficiently, the latent 
Gaussian model has to satisfy some properties. The latent 
Gaussian field x, generally of large dimension, must admit 
conditional independence properties. Thus, the latent 
Gaussian field is a Gaussian Markov random field (GMRF) 
with a sparse precision matrix (inverse of covariance matrix) 
Q (Rue and Held 2005), as the efficiency of INLA relies on 
efficient algorithms for sparse matrices. Due to the use of a 
numerical integration scheme and optimization methods 
in INLA, it needs to integrate over the non-Gaussian hyper- 
parameter space, and therefore, the dimension of non- 
Gaussian hyperparameters 0 cannot be too large, say <14. In 
addition, the likelihood for each observation y, depends on 
the latent Gaussian field only through the linear predictor 
= #(/';) > where g(-) is a known link function and fa = £(y,| 
x, 0), that is, n(y\x, 0) = n{yi\r\ h 0). 

The Al (eqn 1) and AZI (eqn 2) animal model can be 
formulated in the INLA framework with a Gaussian likeli- 
hood yi\r\ i :~ M(r] b a 2 ) and an identity link function, r\ { = 
jU;, where f? ; is the linear predictor. The linear predictor in 
the Al model can be written as: 

rii = po+zfp + Oi, (3) 
and the linear predictor in the AZI model as: 

>h = Po+zlP + a i + z i- W 

It is shown (Henderson 1976; Quaas 1976; Steinsland 
and Jensen 2010; Holand et al. 2013) that the inverse of 
the autosomal relationship matrix A -1 is a sparse matrix, 
which can be calculated from the pedigree. Further, the 
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inverse of the sex-linked relationship matrix Z^ 1 is also a 
sparse matrix, which can be calculated from the pedigree 
and sex information (Fernando and Grossman 1990). The 
autosomal and Z-linked genes are on different chromo- 
somes; therefore, a and z are assumed independent, and 
their joint precision matrix is also sparse. These two pre- 
cision matrices are easily fitted into the INLA framework. 
The latent field x = (/?, a, z) therefore admits conditional 
independence properties, such that x is a GMRF, where 
the precision matrices for the latent field are sparse. As 
the number of non-Gaussian hyperparameters 
0 = {a 2 ^a 2 al a 2 ,a 2 ) is small, and the likelihood of each 
observed trait, depends on the latent field only through 
the linear predictor r\ b the requirements for INLA are ful- 
filled also for the AZI animal model. 

The R software (R Development Core Team 2013) were 
used in our study. The R-INLA package (available at: 
http://www.r-inla.org) makes inference from GRMF mod- 
els using the INLA methodology. Further, the R-package 
AnimallNLA includes functionality for calculating A -1 
and Z~ , and can be downloaded at: http://www.r-inla. 
org. The details of the procedure to efficiently compute 
Z~ directly from pedigree and sex information are given 
in the Appendix SI. 

Model comparison 

Model comparisons in both the simulation study and in 
the barn owl case study are carried out using the deviance 
information criterion (DIC), which is a measure of com- 
plexity and fit (Spiegelhalter et al. 2002). The model with 
the smallest DIC is considered the best model, and 
according to Spiegelhalter et al. (2002), differences in 
DIC, ADIC, of more than 10 should definitely rule out 
the model with the higher DIC. 

In Holand et al. (2013), a simulation-based test of the 
ability of the difference in DIC to chose between animal 
models with and without genetic effects was presented. 
Here, we followed the same ideas and conducted a simu- 
lation-based hypothesis test to test whether Z-linked 
inheritance can be identified using ADIC. Under the null 
hypothesis, H 0 , the Al animal model is true. Under the 
alternative hypothesis, Hi, the AZI animal model is true. 
To estimate the probability of type-I error (reject H 0 
when it is true), we sample S datasets from the Al animal 
model. For each of these s = 1,. . ., S datasets, we fit both 
an Al model and an AZI animal model and calculate the 
difference in DIC, ADIQ = DIC(AI) S -DIC(AZI) S . The 
obtained S values of ADIC are then be used as an approx- 
imation to the sampling distribution of ADIC under the 
null hypothesis. As we reject the null hypothesis for 
ADIC > 10, the proportion of ADIC > 10 is an estimate 
for the probability of type-I error. 
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We also find the power (the probability of rejecting H 0 
when H 1 is true) of the test for some chosen values of a 2 , 
a 2 , and a 2 . For each parameter set, we sample S datasets 
from the AZI animal model, fit both an AI animal model 
and an AZI animal model, and calculate ADIC for these 
models. The proportion of the S ADIC values larger than 
our chosen limit of ADIC =10 is an estimate of power 
when there are some sex-linked genetic effects. 

Simulation study 

The aim of the simulation study was threefold: (1) to 
assess the impact of ignoring Z-linked inheritance on esti- 
mated variance components when sex-linked inheritance 
is present, (2) to evaluate the ability of ADIC to choose 
between models with and without sex-linked genetic 
effects, and (3) to evaluate bias and coverage of variance 
parameters for the two animal models (AI and AZI). 

For the simulation study, we use the barn owl pedigree 
presented in Material and Methods, Field data, and we 
also impose the same missing data structure in the simu- 
lated dataset as in the barn owl dataset. Therefore, we can 
also validate whether the barn owl study system is suitable 
for identifying Z-linked inheritance. We sample data from 
the AZI animal model defined in eqn (2) for the pedigree 
described in Material and Methods, Animal model for 
chosen sets of parameters. These parameters are chosen as 
follows: we simulate approximately standard Gaussian 
datasets by setting ft = 0 and the total variance 
a\ + (3/4) + a] = l. Further, the heritability 



is fixed to h 2 = 0.6, hence a\ = 0.4 and 
h 2 = a 2 a + 3/4(7:;. By choosing <s\ = {0, 0.1, 0.2, 0.3, 0.4, 
0.5, 0.6, 0.7, 0.8}, the corresponding values for autosomal 
variance are a 2 a = {0.6, 0.525, 0.450, 0.375, 0.3, 0.225, 
0.150, 0.075, 0}. These parameter sets range from only 
autosomal inheritance (o^o^cf) = (0.6,0,0.4), that is, 
the AI animal model, to only sex-linked inheritance 
(*»') = (0,0.8,0.4). 

One thousand replicated datasets (S = 1000) were sim- 
ulated for each of these nine parameter sets. Each dataset 
was fitted to both the AI animal model in eqn (1) and 
the AZI animal model in eqn (2). Posterior mean and 
95% credible intervals for variance parameters as well as 
DIC were calculated for each model, and ADIC for each 
pair of models. 

To summarize the simulation results, we also calculated 
the bias and coverage for each parameter set and each 
model. Bias is a measure of the accuracy of an estimator 
9 of a parameter 9, and defined as Bias(fJ) = E(0) — 9, 
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where E(9) = 9 is the average of the 1000 estimated vari- 
ance parameters and 9 is the parameter value used in the 
simulations. Further, we use the coverage to assess the 
precision of the estimator, which is the proportion of 
times the true parameter 9 falls within the 95% credible 
interval of 9, calculated as number of times the parameter 
value used in the simulations are either larger or smaller 
than the estimated lower (0.025%) and upper (0.975%) 
quantiles out of the 1000 simulations. 

Model for field data 

We analyze the data in Materials and Methods, Field data, 
where the plumage spot diameter is assumed to have a 
Gaussian likelihood using the AZI animal models defined 
in Materials and Methods, Animal model. The inference 
is carried out using INLA described in Materials and 
Methods, Modeling Z-linked inheritance - INLA and 
computational issues. 

First, we do a model comparison using DIC to choose 
which fixed effects (sex and hatch year) and random 
effects (autosomal or Z-linked additive genetic effect) to 
include in our model. 

We started with a full model: y { = /? 0 + ftex(i) + thatch year 
(i) + a,- + z,- + and removed one variable at a time in a 
stepwise manner. In each step, all nested models are 
examined, where we only report the one with the lowest 
DIC (i.e., the best at each step). For comparison, we also 
fitted the best model without sex-linked variance. For all 
fixed-effects parameters [3, we use the prior for the covariates 
j3 t ~jV(0, 100). 

To examine whether any evolutionary trends in poster- 
ior mean of mean additive genetic effects in spot diameter 
had occurred over the study period, we find linear combi- 
nations of both autosomal and Z-linked mean additive 
genetic effects for each hatch year (i.e., cohort); 

EfeC ye>r ( 1 / N 'year)fl, and Eiec^CVNyear)*;, where N y™r is 
the number of individuals with a given hatch year and 
summing over all these individuals (Sorensen et al. 1994). 
For further details, see Holand et al. (2013). 

Results 

Simulation study 

To evaluate the ability of ADIC to identify sex-linked 
inheritance, we consider the results of ADIC from the 
simulation study, see Fig. 1 panel A where boxplots of 
obtained ADIC for different values of a\ are plotted. Of 
the S = 1000 datasets simulated from a model with only 
autosomal inheritance {a 2 = 0), we find that only eleven 
datasets have ADIC > 10, and hence, we have an esti- 
mated probability of type-I error (i.e., significance level) 
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of only 0.011. The power of the test (i.e., the probability 
of correctly rejecting H 0 when there is sex-linked inheri- 
tance) can be found in Fig. 1 panel B as a function of a\. 
We find that for g\ = 0.1, the power is only 0.23, but it 
increases fast and is already 0.53 for a\ = 0.2 and 0.73 
for al = 0.3. Hence, for the barn owl system, we are able 



to detect sex-linked inheritance if there is a relatively sub- 
stantial amount of sex-linked effects. 

To evaluate the consequences of not including sex- 
linked inheritance in the model when it is present, we con- 
sider Fig. 1 panel C where the estimated values of a\ are 
plotted when fitting an AI animal model against the true 
value (gray lines). We find that regardless of the true value 
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Figure 1. Results from simulation studies showing performance of animal models for estimating Z-linked additive genetic effects. (A) Boxplot of 
simulated values of ADIC (limit ADIC = 10 indicated as a horizontal solid line) against the value of cr| (AZI) used in the simulations (together with 
(jj and al). (B) Posterior mean (filled squares/solid lines) with 95% credible interval (dashed line) for a\ (AZI) from the simulation study (together 
with a\ and o%), power of the model selection test using ADIC >10 as limit (x'es/solid line) estimated using the simulation approach, and a 1:1 
function of true vs. estimated parameter values (gray line). (C) Posterior mean (open triangles/solid line) for a\ (Al) (gray) with 95% credible 
interval (dashed lines) and for a\ (AZI) (open squares, black) with 95% credible interval (dashed lines) from the simulation study (together with a\ 
and al), power of the model selection test, and 1:1 function as described in panel (A). (D) Posterior mean (solid lines) for al (Al) (gray) with 95% 
credible interval (dashed lines) and for al (AZI) (black) with 95% credible interval (dashed lines) from the simulation study against the value of al 
used in the simulations (together with al) and power function as described in (A). 
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of c\, the estimate is close to the total amount of additive 
variance {g 1 + 3/4<x^ = 0.6). This results in large bias 
and low coverage for a 2 a when fitting an AI animal model 
when sex-linked inheritance is present, see Table SI. From 
Fig. 1 panel D, we see that not including sex-linked inheri- 
tance has very little effect on the estimated values of a 2 . 

From Table SI we see that, when fitting the AZI animal 
model, the bias is small and coverage is good except for 
small values of the additive variances a 2 a and a 2 . This is 
known to be due to prior sensitivity (see Holand ef al. 
2013). When models are fitted to a dataset without sex- 
linked inheritance {a\ = 0 and a 2 = 0.6 in out simula- 
tion study), we see from Fig. 1 panel B, C, and D that the 
AZI estimates perform slightly worse than the AI esti- 
mates in terms of larger credible intervals. 

Plumage spot diameter in barn owls 

The results of the model comparison where different models 
were fitted to the spot diameter data and compared in terms 
of their DIC values are fisted in Table 1. The difference 
in DIC values for the model, including both autosomal and 
Z-linked effects versus the model accounting only for 
autosomal effects, was 180 in favor of the model which 
explicitly models Z-linked additive genetic effects. Difference 
in DIC thus greatly exceeds the chosen limit of 10 and 
decisively indicates that spot diameter is partially influenced 
by variation in Z-linked genes. Furthermore, the overall 
best model includes sex as fixed effect in addition to both 
autosomal and Z-linked additive genetic effects. 

The estimated posterior mean and 95% credible inter- 
val for the additive genetic variances were a\ = 0.4405 
(95% CI: 0.3603 to 0.5336), a\ = 0.2719 (95% CI: 
0.1833 to 0.3880), and a] = 0.2012 (95% CI: 0.1639 to 
0.2439). The parameter a 2 is the Z-linked additive genetic 
variance for noninbred males. Thus, for females, we have 
& = 0.1360 (95% CI: 0.0917 to 0.1940). 

In comparison, a model including sex and only autoso- 
mal additive genetic effect yields additive genetic variance 
a\ = 0.6124 (95% CI: 0.5392 to 0.6936) and 
al = 0.2199 (95% CI: 0.1799 to 0.2650). 



Table 1. Different model specifications explaining variance in spot 
diameter of Swiss barn owls and the corresponding deviance informa- 
tion criteria (DIC). 



Model 



DIC 



Sex + hatch year + autosomal effect + Z-linked effect 5050 

Sex + autosomal effect + Z-linked effect 4757 

Sex + autosomal effect 4937 

Autosomal effect 5898 



The linear combinations of posterior mean of mean 
additive genetic effects across cohorts suggest for spot 
diameter that there was an increase in additive genetic 
effect for autosomal loci, but no increase in the additive 
genetic effect for Z-linked loci (Fig. 2). To test this, we 
investigated whether the difference between cohorts 
1996 and 2007 in posterior mean additive genetic effects 
was significant for either autosomal and Z-linked loci: 
J2i€C m6 (1/Ni9 96 )fli - E;ec 2M7 (1/ N 2007)«; and £i 6 c 1996 



(A) 



autosomal 
z-linked 




2000 2002 

Hatch year 



(B) 



£ 

w 
o 
Q_ 




The best model is given in bold. 



Figure 2. (A) Posterior mean of mean additive genetic effect for spot 
diameter of all individuals in the pedigree for each cohort (i.e., hatch 
year) 1996-2007 for autosomal loci (black) and Z-linked loci (gray) 
(solid lines) with 95% credible intervals (dashed lines). The mean spot 
diameter was standardized to have mean 0 and variance 1 . (B) 
Posterior of difference between cohorts 1996 and 2007 in mean 
additive genetic effects for autosomal loci (solid lines) and Z-linked 
loci (dashed lines) for spot diameter in Swiss barn owls. 
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(1/N 1996 )z; - J2i € c 2mj (V N 2007)z;, respectively. The dif- 
ference between additive genetic effects for cohorts 1996 
and 2007 was significant for autosomal loci, with mean 
difference -0.206 (SD = 0.055) and 95% CI 
(—0.313,-0.097). In contrast, the difference was not sig- 
nificant for Z-linked loci; mean difference —0.052 
(SD = 0.050) and 95% CI (-0.152,0.046). The posterior 
marginals of the difference for autosomal and Z-linked 
loci are given in Fig. 2. 

Discussion 

From the simulation results, we see that although we have 
modeled only autosomal inheritance, in the presence of 
Z-linked inheritance in a trait, the total amount of addi- 
tive genetic variance is correct. This is apparently because 
all the additive genetic variance, including the part due to 
genes on the Z-chromosome, is picked up by the esti- 
mated autosomal additive genetic variance. Hence, using 
an AI animal model when an AIZ model is true gives an 
estimate of autosomal additive genetic variance, which 
corresponds to the total amount of additive genetic effects 
in the AZI model. Incorrect predictions of responses to 
selection may, however, be one of the consequences of 
not modeling Z-linked inheritance when it exists. Results 
from the simulation study therefore clearly illustrate the 
importance of specifying the correct model in the pres- 
ence of Z-linked inheritance. Further, the simulation 
study also demonstrates that for our barn owl study sys- 
tem (i.e., this type of pedigree and missing data struc- 
ture), difference in DIC between AI and AIZ models is a 
good measure for model choice. Study systems that show 
low support for sex-linked effects, for example, Husby 
et al. (2013), would benefit from a simulation study to 
explore the system's ability to separate autosomal and 
sex-linked additive genetic effects. In our simulation 
study, we used 1000 datasets; however, more datasets 
might be needed if a more robust estimation of a specific 
power is desired. 

The analysis of spot diameter in the empirical barn 
owl dataset showed that spot diameter is clearly influ- 
enced by both autosomal and Z-linked additive genetic 
effects. The results show that spot diameter is under 
strong genetic influence, and genes on the Z-chromo- 
some contribute a substantial amount to the total phe- 
notypic variation. 

According to theory on effects of selection when it 
mainly acts on the heterogametic sex, for example, the 
females in birds and males in mammals, we expect to see 
a change in mean additive genetic effects mainly in the 
genes found in the autosomes (Charlesworth et al. 1987). 
This is in accordance with the results found in this study, 
where the changes in posterior mean of mean additive 



genetic effects of spot diameter across cohorts suggest an 
increase in autosomal additive genetic effects over the 
study period (Fig. 2A). This is supported by the signifi- 
cant difference between cohort 1996 and 2007 found for 
the autosomal additive genetic effects, whereas there was 
no significant change across cohorts for the Z-linked 
additive genetic effects (Fig. 2B). However, it is difficult 
to determine whether the observed change in mean 
breeding values is due to an evolutionary response to 
selection on spot diameter or random genetic drift, as 
genetic drift may cause independent fluctuations in breed- 
ing values across generations (Hadfield et al. 2010). In 
any case, the result that the genetic changes mainly 
occurred on the autosomes corresponds with other stud- 
ies of birds, suggesting that most of sexually antagonistic 
genes beneficial for females are located on the autosomes 
(Ellegren and Parsch 2007; Mank and Ellegren 2009). 

Another possible explanation to the small change in 
breeding values over the cohorts for Z-linked genes is that 
the spot diameter itself is not under selection, but rather 
another trait that is genetically correlated with spot dia- 
meter on the autosomal chromosomes is under selection. 
This is in accordance with findings in Roulin and Ducrest 
(2011), which showed that spot size displayed by mothers 
is correlated with offspring quality measures including 
parasite resistance, resistance to oxidative stress, and an 
increase in corticosterone levels, appetite, and the ability 
to withstand lack of food. 

Sex-specific selection is the process in which selection 
is favoring different optimal character states in the two 
sexes, a mode of selection that recently has received much 
attention by evolutionary biologists (see e.g., Lande 1980; 
Foerster et al. 2007; Cox and Calsbeek 2009; Mills et al. 
2012; Stearns et al. 2012). Modeling sex-linked genetic 
variance and performing a simulation study to explore 
the strength of the study system to identify sex-linked 
genetic variance are especially important when working 
with sexual conflict and sex-specific selection. The covari- 
ance between a given trait and selection can be positive 
for males and negative for females, or the other way 
around. This type of selection may, for example, occur 
because the two sexes have differing roles in reproduc- 
tion, leading to different phenotypic optima in females 
and males. The study of sex-specific selection is interest- 
ing because this pattern of selection may account for the 
evolutionary stability of sexual dimorphism, it may also 
explain why genetic variation is not eroded, and it pro- 
vides interesting implications into the understanding of 
intralocus genetic conflict (Bonduriansky and Chenoweth 
2009). This type of conflict results from the fact that dif- 
ferent alleles are favored in the two sexes, which can 
result in intricate phenomena such as sex ratio bias 
(Blackburn et al. 2010). 



2014 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd. 



1563 



Sex-Linked Inheritance and Bayesian Animal Model 



C. T. Larsen ef al. 



Our models assume that the additive genetic effects are 
the same in both sexes, that is, the intersex genetic corre- 
lation is one. As the spot diameter in the Swiss barn owls 
is a sexually dimorphic trait and has been shown to be 
subject to sex-specific selection (Roulin et al. 2010), this 
assumption could be violated, which might results in 
some bias in the estimated autosomal and sex-linked 
additive genetic variances. To explore this, a model treat- 
ing the spot diameter in males and females as two differ- 
ent traits that are genetically correlated has to be fitted. 
This is, however, outside the scope of this study, but 
studies of different species suggest that intersex correla- 
tions of homologous morphological traits often are close 
to one (Jensen et al. 2008; Kruuk et al. 2008). 

The animal model we have used (implicitly) assumes 
that both autosomal and sex-linked genetic effects are 
additive, that there is no sex chromosome dosage com- 
pensation, and that missing observations are missing at 
random. Nonadditive effects such as dominance and epis- 
tasis are known to be hard to identify from nonexperi- 
mental study systems (Lynch and Walsh 1998), and it is 
outside the scope of this work to our extend model to 
account for these effects. Other studies have suggested 
that sexually antagonistic genes often are dominant (Elle- 
gren and Parsch 2007; Mank and Ellegren 2009). How- 
ever, the importance of nonadditive genetic effects is 
arguable as some studies suggest both that dominance 
and epistasis effects may contribute little to the pheno- 
typic variance (Merilii et al. 2001; Visscher et al. 2007; 
Crow 2010, but see e.g. Carlborg and Haley 2004). 

Hence, dominance and epistasis may not affect our 
results considerably. The assumption of no sex chromo- 
some dosage compensation seems to be a good assump- 
tion as an overall dosage compensation has not been 
found in birds (Ellegren et al. 2007a; Itoh et al. 2007). 
The assumption that missing observations of the trait of 
interest are missing at random is further explored in 
Steinsland et al. (2014), where it is concluded that for this 
system, the assumption does not influence the variance 
estimates to any large extent. 

Both the model and the data used in Roulin et al. 
(2010) and our study are slightly different. Therefore, we 
do not expect results in these two papers to coincide 
exactly. However, the results found in the two studies are 
essentially identical when it comes to sex-linked variances. 
In Roulin et al. (2010), the Z-linked additive genetic vari- 
ance a\ was reported for females (a\ = 0.132), while it is 
reported for males here. Consequently, if we compare 
(2 x olj = 0.264) with a\ = 0.27, we see that the addi- 
tive sex-linked variance is very similar. Furthermore, the 
trends in additive genetic effects are similar in Roulin 
et al. (2010) and in the current paper. Posterior distribu- 
tions of mean difference in mean breeding values were, 



however, not exactly equal in the two studies, probably as 
a consequence of different models and data used. In Rou- 
lin et al. (2010), phenotypes of more owls were included, 
data were standardized within each sex (see also Steins- 
land et al. 2014), and hatch year but not sex was included 
as fixed effect in their animal model. It should be noted 
that standardizing the phenotypes within each sex forces 
the variance within each sex to be equal, while an animal 
model with sex-linked effects implicitly assumes that 
males have larger variance than females. Hence, it is 
inconsistent to do sex-specific standardization prior to 
applying an animal model with sex-linked inheritance. 
Finally, different methods for computing Z~ l were used. 
In Roulin et al. (2010), the software Mendel (Lange et al. 
2001) was used to compute Z, and MATLAB to invert 
this matrix to obtain Z~ x . Numerical problems with this 
approach were reported. 

Conclusion 

We have in this study introduced a methodology for esti- 
mation and testing identifiability issues regarding sex- 
linked additive genetic effects and discussed consequences 
of not modeling this variance when it is present. Through a 
simulation study, we have shown that for a real wild popu- 
lation system (with a given pedigree, missing data structure, 
and sex distribution) that both autosomal and sex-linked 
effects can be estimated, these effects can be distinguished 
(i.e., they are identifiable), and difference in DIC between 
animal models with only autosomal inheritance and both 
autosomal and sex-linked inheritance can be used to test 
whether sex-linked inheritance is present. Using an animal 
model with only autosomal inheritance when sex-linked 
inheritance is present results in inflated estimates of the 
autosomal additive variance, as it also includes the sex- 
linked variance of the population. This might give mislead- 
ing interpretations, especially when response to sex-specific 
selection is studied, as the heterogametic sex for instance 
will have a slower response to selection than the homoga- 
metic sex for genes on the sex chromosome when alleles 
have largely dominant effects (Charlesworth et al. 1987). 
We are not able to obtain any knowledge about potential 
sex-linked inheritance from a model assuming only autoso- 
mal inheritance. On the other hand, fitting an animal 
model with both autosomal and sex-linked effects to a sys- 
tem where no sex-linked effects are present, performs 
approximately equally well as the model that (correctly) 
assumes only autosomal inheritance. We therefore recom- 
mend that animal models including both autosomal and 
sex-linked effects are used, or at least tested. 

In our study of plumage spot diameter in a Swiss barn 
owl population, we found that sex-linked effects account 
for a substantial proportion of the phenotypic variance. 
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Earlier results indicated that this trait is under sex-specific 
selection (Roulin et al. 2010). 
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